/* 
 * File:   LocalSurfaceCalculator.h
 * Author: boyan
 *
 * Created on August 26, 2009, 1:46 PM
 */

#ifndef _LOCALSURFACECALCULATOR_H
#define	_LOCALSURFACECALCULATOR_H

#include <vector>

#include <Vector3D.hpp>
using Mesquite::Vector3D;

#include <Matrix3D.hpp>
using Mesquite::Matrix3D;

// **************************************************************************

class LocalSurfaceCalculator
{
protected:
    Vector3D FiducialPoint;
    Vector3D t, s, n; // the local basis;

public:

    LocalSurfaceCalculator() :
    FiducialPoint(),
    t(1, 0, 0),
    s(0, 1, 0),
    n(0, 0, 1)
    {
    };

    LocalSurfaceCalculator(const LocalSurfaceCalculator &to_copy) :
    FiducialPoint(to_copy.FiducialPoint),
    t(to_copy.t),
    s(to_copy.s),
    n(to_copy.n)
    {
    };

    LocalSurfaceCalculator & operator=(const LocalSurfaceCalculator &to_copy)
    {
        FiducialPoint = to_copy.FiducialPoint;
        t = to_copy.t;
        s = to_copy.s;
        n = to_copy.n;
        return *this;
    }

    Vector3D &getFiducialPoint();
    const Vector3D &getFiducialPoint() const;
    void setFiducialPoint(const Vector3D &pt);

    void setLocalBasisVectors(const Vector3D &t, const Vector3D &s, const Vector3D &n);
    void setLocalBasisXYZ();
    Vector3D &getLocalBasisByIndex(int index);
    const Vector3D &getLocalBasisByIndex(int index) const;

    virtual void flipNormal();

    virtual void computeLocalSurface() = 0;
    virtual Vector3D getNormal() const = 0; // at FiducialPoint
    virtual Vector3D getNormal(const Vector3D &pt) const = 0; // at projection of pt
    virtual Vector3D getProjection(const Vector3D &pt) const = 0;

};

// **************************************************************************

class PointSet : public std::vector< Vector3D >
{
public:
    Vector3D MeanPoint() const;
    Matrix3D ScatterMatrix() const;
    Matrix3D ScatterMatrix(const Vector3D & pt) const;
    void ComputeLocalBasis(Vector3D &t, Vector3D &s, Vector3D & n) const;
    void ComputeLocalBasisAtPoint(const Vector3D &pt, Vector3D &t, Vector3D &s, Vector3D & n) const;
    bool FitSphereAtPoint(const Vector3D &pt, Vector3D &centre, double &radius) const;
    void FitPlane(Vector3D &point, Vector3D &normal) const;
    void FitPlaneAtPoint(const Vector3D &point, Vector3D &normal) const;
};

class PointSetSurfaceCalculator : public LocalSurfaceCalculator
{
public:
    PointSet AdjacentPoints;

    PointSetSurfaceCalculator() :
    LocalSurfaceCalculator(),
    AdjacentPoints()
    {
    };

    PointSetSurfaceCalculator(const PointSetSurfaceCalculator &to_copy) :
    LocalSurfaceCalculator(to_copy),
    AdjacentPoints(to_copy.AdjacentPoints)
    {
    };

    PointSetSurfaceCalculator & operator=(const PointSetSurfaceCalculator &to_copy)
    {
        LocalSurfaceCalculator::operator =(to_copy);
        AdjacentPoints = to_copy.AdjacentPoints;
        return *this;
    }

    // used to allocate copies of derived classes
    virtual PointSetSurfaceCalculator *NewCopy() const = 0;
#define implement_NewCopy(TYPE)   \
    virtual PointSetSurfaceCalculator *NewCopy() const {return new TYPE(*this);}

};

class LocalBall : public PointSetSurfaceCalculator
{
protected:
    double Kappa;
    Vector3D BallCentre;
    signed int NormalSign;

public:

    LocalBall() :
    PointSetSurfaceCalculator(), // this one is called by default
    Kappa(0.0),
    BallCentre(), // default constructor sets it to (0,0,0)
    NormalSign(1)
    {
    };

    LocalBall(const LocalBall &to_copy) :
    PointSetSurfaceCalculator(to_copy),
    Kappa(to_copy.Kappa),
    BallCentre(to_copy.BallCentre),
    NormalSign(to_copy.NormalSign)
    {
    };

    LocalBall & operator=(const LocalBall &to_copy)
    {
        PointSetSurfaceCalculator::operator =(to_copy);
        Kappa = to_copy.Kappa;
        BallCentre = to_copy.BallCentre;
        NormalSign = to_copy.NormalSign;
        return *this;
    }

    implement_NewCopy(LocalBall);

    virtual void flipNormal()
    {
        if(NormalSign)
            NormalSign *= -1;
        else
            NormalSign = -1;
    }

    virtual void computeLocalSurface();
    virtual Vector3D getNormal() const; // at FiducialPoint
    virtual Vector3D getNormal(const Vector3D &pt) const; // at projection of pt
    virtual Vector3D getProjection(const Vector3D &pt) const;
};

// **************************************************************************

inline Vector3D &LocalSurfaceCalculator::getLocalBasisByIndex(int index)
{
    static Vector3D nilvec;
    switch(index)
    {
        case 0: return this->t;
        case 1: return this->s;
        case 2: return this->n;
        default: return nilvec;
    }
}

inline const Vector3D &LocalSurfaceCalculator::getLocalBasisByIndex(int index) const
{
    static Vector3D nilvec;
    switch(index)
    {
        case 0: return this->t;
        case 1: return this->s;
        case 2: return this->n;
        default: return nilvec;
    }
}

inline Vector3D &LocalSurfaceCalculator::getFiducialPoint()
{
    return this->FiducialPoint;
}

inline const Vector3D &LocalSurfaceCalculator::getFiducialPoint() const
{
    return this->FiducialPoint;
}

inline void LocalSurfaceCalculator::setFiducialPoint(const Vector3D &pt)
{
    this->FiducialPoint = pt;
}

inline void LocalSurfaceCalculator::setLocalBasisXYZ()
{
    this->t = Vector3D(1, 0, 0);
    this->s = Vector3D(0, 1, 0);
    this->n = Vector3D(0, 0, 1);
}

inline void LocalSurfaceCalculator::setLocalBasisVectors(const Vector3D &t,
        const Vector3D &s, const Vector3D &n)
{
    this->t = t;
    this->s = s;
    this->n = n;
}

inline void LocalSurfaceCalculator::flipNormal()
{
    this->n *= -1;
    this->s = this->n * this->t; // keep positive orientation
}

// ***************************************************************************

inline Matrix3D PointSet::ScatterMatrix() const
{
    return ScatterMatrix(MeanPoint());
}

inline void PointSet::ComputeLocalBasis(Vector3D &t, Vector3D &s, Vector3D &n) const
{
    ComputeLocalBasisAtPoint(MeanPoint(), t, s, n);
}

#endif	/* _LOCALSURFACECALCULATOR_H */

